import matplotlib.pyplot as plt
import numpy as np
#函数1
def func1(x,y):
    return (5*y-6*np.exp(-x))
#问题1的准确解
def func(x):
    return (np.exp(-6*x)+1.0e-13)*np.exp(5*x)
#欧拉法
def ola(x,y0,h):
    y=np.zeros(len(x))
    #初值
    y[0]=y0
    for i in range(len(x)-1):
        #根据上一个值求下一个值
        y[i+1]=y[i]+h*func1(x[i],y[i])
        
    return y
#三阶龙格-库塔法
def lunge_3(x,y0,h):
    y=np.zeros(len(x))
    y[0]=y0
    k1,k2,k3=0,0,0
    for i in range(len(x)-1):
        k1=func1(x[i],y[i])
        k2=func1(x[i]+h/2,y[i]+h*k1/2)
        k3=func1(x[i]+h,y[i]-h*k1+2*h*k2)
        y[i+1]=y[i]+(k1+4*k2+k3)*h/6

    return y
#四阶龙格-库塔法
def lunge_4(x,y0,h):
    y=np.zeros(len(x))
    y[0]=y0
    k1,k2,k3,k4=0,0,0,0
    for i in range(len(x)-1):
        k1=func1(x[i],y[i])
        k2=func1(x[i]+h/2,y[i]+h*k1/2)
        k3=func1(x[i]+h/2,y[i]+h*k2/2)
        k4=func1(x[i]+h,y[i]+h*k3)
        y[i+1]=y[i]+(k1+2*k2+2*k3+k4)*h/6
        
    return y
#问题1
fx0=1.0+1.0e-13
x1=np.arange(0,6,0.0001)
y1=ola(x1,fx0,0.0001)
y2=lunge_3(x1,fx0,0.0001)
y3=lunge_4(x1,fx0,0.0001)
fig,ax=plt.subplots(2,2,figsize=(18,6))
y_=func(x1)
loss1=y_-y1
loss2=y_-y2
loss3=y_-y3
print("欧拉法各点的误差为：\n",loss1)
print("欧拉法各点的误差为：\n",loss2)
print("欧拉法各点的误差为：\n",loss3)
print("平均误差为分别为：\n","ola:",sum(loss1)/len(loss1),"lunge3:",sum(loss2)/len(loss2),"lunge4:",sum(loss3)/len(loss3))
ax[0][0].plot(x1,y1,c='r',label="数值解")
ax[0][1].plot(x1,y2,c='r',label="数值解")
ax[1][0].plot(x1,y3,c='r',label="数值解")
ax[1][0].legend()
ax[1][1].plot(x1,y_,c='b',label="精确解")
ax[1][1].legend()
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus'] = False
plt.show()